/*
 * HaliteSolubility.cpp
 *
 *  Created on: 12 Sep 2011
 *      Author: Allan
 */

#include "HaliteSolubility.h"

// C++ includes
#include <cmath>
#include <iostream>
#include <iomanip>
using namespace std;

// GeoReact includes
#include "../Assembler.h"
#include "../Applications/Equilibrator.h"
#include "../Utilities/Range.h"

const double HaliteMoles = 10.0;

double HaliteSaturation(double T, double P);

void DemoHaliteSolubility()
{
	System system = AssembleSystem({"H2O", "Halite"});
	
	Equilibrator equilibrator;
	
	cout << setw(15) << left << "T";
	cout << setw(15) << left << "Solubility";
	cout << setw(15) << left << "Na[+]";
	cout << setw(15) << left << "Sawamura";
	cout << setw(15) << left << "gNa[+]";
	cout << setw(15) << left << "gCl[-]";
	cout << setw(15) << left << "gNaCl(aq)";
	
	cout << endl;

	BOOST_FOREACH(double T, Range(25.0, 100.0, 1.0))
	{
		system.SetTemperature(T);
		
		equilibrator.ResetConditions();
		
		equilibrator.AddCondition("ChargeBalanced");
		equilibrator.AddCondition("FixedMass", "H2O", 55.508);
		equilibrator.AddCondition("FixedTotalMass", "Na[+]", HaliteMoles);
		equilibrator.AddCondition("FixedTotalMass", "Cl[-]", HaliteMoles);
		
		equilibrator.Equilibrate(system);
		
		cout << setw(15) << left << T;
		cout << setw(15) << left << HaliteMoles - system.Get("Moles", "Halite");
		cout << setw(15) << left << system.Get("Moles", "Na[+]");
		cout << setw(15) << left << HaliteSaturation(T, 1.0);
		cout << setw(15) << left << system.Get("ActivityCoefficient", "Na[+]");
		cout << setw(15) << left << system.Get("ActivityCoefficient", "Cl[-]");
		cout << setw(15) << left << system.Get("ActivityCoefficient", "NaCl(aq)");
		cout << endl;
	}
}

double HaliteSaturation(double T, double P)
{
	double Tk = T + 273.15;
	
	return exp(2.86623 - 571.361/Tk + 77128/Tk/Tk + 
			   4.0479E-5*P - 7.445E-9*P*P + 6.209E-13*P*P*P);
}
